#this program demonstrates some of the network randomization bits and 
#example models.

#clear everything to start...
rm(list = ls())
gc()
##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 375352 20.1     750400 40.1   592000 31.7
## Vcells 583261  4.5    1308461 10.0   786378  6.0
#load basic data manipulation bits
library(dplyr); 
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(readr);
library(magrittr)     #Supports pipe (%>%) commands that allow you to perform multiple operations with one statement
library(tidyr)        #Additional functions for manipulating data
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
## 
##     extract
library(ggplot2)

#load the data - this is the same code as before
setwd("/Users/mariacristinaramos/Documents/MARIA CRISTINA/ACTIVE/SN-H2018/")
AHS_Base <- read_csv ("ahs_wpvar (1).csv",
                      col_names = TRUE);
## Parsed with column specification:
## cols(
##   .default = col_integer(),
##   POP_BEH = col_double(),
##   PSICK_SS = col_double(),
##   PSICK_H = col_double(),
##   PSICK_L = col_double(),
##   POP_BEH_SS = col_double()
## )
## See spec(...) for full column specifications.
AHS_adjlist <- AHS_Base %>%
  select(ego_nid, mfnid_1:mfnid_5, ffnid_1:ffnid_5, grade, sex, PSMOKES,commcnt) %>%
  filter(commcnt==1);

AHS_Edges <- AHS_adjlist %>%
  rename( id = `ego_nid`,
          gender = `sex`) %>%
  gather(Alter_Label, Target, mfnid_1:mfnid_5, ffnid_1:ffnid_5, na.rm = TRUE)

AHS_Edges=AHS_Edges %>% filter (Target != 99999);
AHS_Edges=AHS_Edges %>%select(id, Target);

#install.packages("statnet")
library(statnet)
## Loading required package: tergm
## Loading required package: statnet.common
## 
## Attaching package: 'statnet.common'
## The following object is masked from 'package:base':
## 
##     order
## Loading required package: ergm
## Loading required package: network
## network: Classes for Relational Data
## Version 1.13.0 created on 2015-08-31.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##                     Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
##                     Martina Morris, University of Washington
##                     Skye Bender-deMoll, University of Washington
##  For citation information, type citation("network").
##  Type help("network-package") to get started.
## 
## ergm: version 3.8.0, created on 2017-08-18
## Copyright (c) 2017, Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
##                     Carter T. Butts, University of California -- Irvine
##                     Steven M. Goodreau, University of Washington
##                     Pavel N. Krivitsky, University of Wollongong
##                     Martina Morris, University of Washington
##                     with contributions from
##                     Li Wang
##                     Kirk Li, University of Washington
##                     Skye Bender-deMoll, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").
## NOTE: Versions before 3.6.1 had a bug in the implementation of the
## bd() constriant which distorted the sampled distribution somewhat.
## In addition, Sampson's Monks datasets had mislabeled vertices. See
## the NEWS and the documentation for more details.
## Loading required package: networkDynamic
## 
## networkDynamic: version 0.9.0, created on 2016-01-12
## Copyright (c) 2016, Carter T. Butts, University of California -- Irvine
##                     Ayn Leslie-Cook, University of Washington
##                     Pavel N. Krivitsky, University of Wollongong
##                     Skye Bender-deMoll, University of Washington
##                     with contributions from
##                     Zack Almquist, University of California -- Irvine
##                     David R. Hunter, Penn State University
##                     Li Wang
##                     Kirk Li, University of Washington
##                     Steven M. Goodreau, University of Washington
##                     Jeffrey Horner
##                     Martina Morris, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("networkDynamic").
## 
## tergm: version 3.4.1, created on 2017-09-12
## Copyright (c) 2017, Pavel N. Krivitsky, University of Wollongong
##                     Mark S. Handcock, University of California -- Los Angeles
##                     with contributions from
##                     David R. Hunter, Penn State University
##                     Steven M. Goodreau, University of Washington
##                     Martina Morris, University of Washington
##                     Nicole Bohme Carnegie, New York University
##                     Carter T. Butts, University of California -- Irvine
##                     Ayn Leslie-Cook, University of Washington
##                     Skye Bender-deMoll
##                     Li Wang
##                     Kirk Li, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("tergm").
## Loading required package: ergm.count
## 
## ergm.count: version 3.2.2, created on 2016-03-29
## Copyright (c) 2016, Pavel N. Krivitsky, University of Wollongong
##                     with contributions from
##                     Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm.count").
## NOTE: The form of the term 'CMP' has been changed in version 3.2
## of 'ergm.count'. See the news or help('CMP') for more information.
## Loading required package: sna
## sna: Tools for Social Network Analysis
## Version 2.4 created on 2016-07-23.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##  For citation information, type citation("sna").
##  Type help(package="sna") to get started.
## 
## statnet: version 2016.9, created on 2016-08-29
## Copyright (c) 2016, Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
##                     Carter T. Butts, University of California -- Irvine
##                     Steven M. Goodreau, University of Washington
##                     Pavel N. Krivitsky, University of Wollongong
##                     Skye Bender-deMoll
##                     Martina Morris, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("statnet").
## unable to reach CRAN
g=as.network(AHS_Edges)
g %v% "grade" <- AHS_adjlist$grade
g %v% "sex" <- AHS_adjlist$sex
g %v% "smokes" <- AHS_adjlist$PSMOKES

# get degree scores
outdegree <- degree(g, cmode = "outdegree")
indegree<- degree(g, cmode = "indegree")
g %v% "indegree" <- indegree
g %v% "outdegree" <- outdegree
g %v% "degree" <- degree(g)
plot(g, vertex.col="grade")

plot(g, vertex.col="sex")

#how transitive is our network?
obs_tran=gtrans(g)
obs_tran
## [1] 0.375
#distribution of dyads
dyads=dyad.census(g)
dyads
##      Mut Asym Null
## [1,]  85  135 2265
#to simulate...first initialize an empty vector
simsize=500;
ran_tran=vector(length=simsize)

#loop for the simulation...
j=1
while(j<(simsize+1)){
  r=rguman(1,nv=71,mut=dyads[1],asym=dyads[2],null=dyads[3])
  ran_tran[j]=gtrans(r)
  if (is.na(ran_tran[j])) next else j=j+1  #skip bad sim cases
}

dt=as.data.frame(ran_tran)
ggplot(dt,aes(x=ran_tran))+geom_histogram()+
  geom_vline(xintercept=obs_tran,
             linetype=1,
             color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#perhaps the clustering observed is due to general categorical mixing, 
#by grade or sex?
grades=as.character(g %v% "grade")
grdmix=mixingmatrix(g, "grade")
grdmix
##        To
## From     7  8  9 10 11 12 Total
##   7     52  5  1  1  0  0    59
##   8      8 33  9  0  1  1    52
##   9      0 10 70  1  4  1    86
##   10     0  0  3 30 10  0    43
##   11     1  0  2  7 43  4    57
##   12     0  0  1  0  2  5     8
##   Total 61 48 86 39 60 11   305
#should be able to do this with the RGNMIX function, but that's 
#not working, so let's do the same thing by hand. 
#here's my failed code for the build in function, 
################
#to simulate...first initialize an empty vector
#simsize=500;
#ran_tran=vector(length=simsize)

#loop for the simulation...
#j=1
#while(j<(simsize+1)){
#  ran_g<-rgnmix(1,grades,grdmix$matrix,method="exact",return.as.edgelist=TRUE)
#  ran_tran[j]=gtrans(as.network(ran_g))
#  j=j+1
#  if (is.na(ran_tran[j])) next else j=j+1  #skip bad sim cases
#}

#dt=as.data.frame(ran_tran)
#ggplot(dt,aes(x=ran_tran))+geom_histogram()+
#  geom_vline(xintercept=obs_tran,
#             linetype=1,
#             color="red")
#now do it by hand.  This code steals directly from some stuff
#jake fisher put together in 2017 session.



#1. get a probability matrix from the mixing matrix...need a denominator
c=as.matrix(table(grades))
cells=c%*%t(c)
mixprob=grdmix$matrix/cells
mixprob
##     To
## From           7           8           9          10          11
##   7  0.520000000 0.038461538 0.025000000 0.006666667 0.000000000
##   8  0.061538462 0.195266272 0.173076923 0.000000000 0.005917160
##   9  0.000000000 0.192307692 4.375000000 0.016666667 0.076923077
##   10 0.000000000 0.000000000 0.050000000 0.133333333 0.051282051
##   11 0.007692308 0.000000000 0.038461538 0.035897436 0.254437870
##   12 0.000000000 0.000000000 0.015625000 0.000000000 0.009615385
##     To
## From          12
##   7  0.000000000
##   8  0.004807692
##   9  0.015625000
##   10 0.000000000
##   11 0.019230769
##   12 0.019531250
#note this is technically not quite right...as we should subtractk
#1 from the diagonal...but I'm being lazy...

# 2. Now create an empty matrix:
prob.matrix <- matrix(NA, ncol = network.size(g),
                      nrow = network.size(g))

# 3. ...fill in the values based on the mixing matrix.  Going to 
#so get the categories...
grade.values <- sort(unique(grades))

# The general idea is to do this...
#prob.matrix[which(grades == 7), which(grades == 7)] <- tie.probs[1, 1]

# ...but that's going to require 36 lines of code.  Let's write a loop instead.
for (i in 1:nrow(mixprob)) {
  for (j in 1:ncol(mixprob)) {
    prob.matrix[which(grades == grade.values[i]), 
                which(grades == grade.values[j])] <- mixprob[i, j]
  }
}

#now we have a matrix where each ij cell is the probabilty of a 
#tie based on mixing...so just do random bernulli draws from that, 
#using the same code strucure as above...

simsize=500;
ran_tran=vector(length=simsize)

#loop for the simulation...
j=1
while(j<(simsize+1)){
  r=rgraph(n = network.size(g), tprob = prob.matrix)
  ran_tran[j]=gtrans(r)
  if (is.na(ran_tran[j])) next else j=j+1  #skip bad sim cases
}

dt=as.data.frame(ran_tran)
ggplot(dt,aes(x=ran_tran))+geom_histogram()+
  geom_vline(xintercept=obs_tran,
             linetype=1,
             color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#could do the same thing by Sex, but the plot we had suggested it's 
#not strongly clustered by sex...but you get the idea...
#also, would be more efficient to create a sample of graphs, store them
#then you can test a bunch of stats agains the same set...but transitivyt
#is cheap to calculate...so using that...

plot(hist(ran_tran),xlim = c(0, .5))

abline(v = obs_tran, col = "red")

# Now we pass both of those to the function in igraph
#one quick test...
 r=igraph::sample_degseq(out.deg = outdegree, in.deg = indegree,method = "simple.no.multiple")
 igraph::transitivity(r)
## [1] 0.1429688
 igraph::plot.igraph(r)

 #same loop bit...could just write a function for this...
simsize=50;
 ran_tran=vector(length=simsize)
 
 #loop for the simulation...note I use the igraph transitivity function
 j=1
 while(j<(simsize+1)){
   r=igraph::sample_degseq(out.deg = outdegree, in.deg = indegree,method = "simple.no.multiple")
   ran_tran[j]=igraph::transitivity(r)
   if (is.na(ran_tran[j])) next else j=j+1  #skip bad sim cases
 }
 
 dt=as.data.frame(ran_tran)
 ggplot(dt,aes(x=ran_tran))+geom_histogram()+
   geom_vline(xintercept=obs_tran,
              linetype=1,
              color="red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

 #again...not even close.  But, all this is binary; let's build a model
 #to look at this with multiple conditionals.
 
 ####################
 # ERGM 
 ####################
 
 #ergm code adapted from Jeff Smith's tutorial.
 #clean house...
 rm(c,cells, dt, dyads, grdmix,prob.matrix,r)
 rm(grade.values,grades,i,j,mixprob,obs_tran,rant_tran,simsize)
## Warning in rm(grade.values, grades, i, j, mixprob, obs_tran, rant_tran, :
## object 'rant_tran' not found
 #super simple model.. 
 mod.rand<-ergm(g~edges)  
## Evaluating log-likelihood at the estimate.
## 
summary(mod.rand)  
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ edges
## 
## Iterations:  6 out of 20 
## 
## Monte Carlo MLE Results:
##       Estimate Std. Error MCMC % p-value    
## edges  -2.7275     0.0591      0  <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 6890  on 4970  degrees of freedom
##  Residual Deviance: 2293  on 4969  degrees of freedom
##  
## AIC: 2295    BIC: 2302    (Smaller is better.)
 #interpretation of model coefficients: 
  #the log-odds of any tie existing is -2.7275
 #probability: exp(-2.7275)/(1+exp(-2.7275))=.0613
exp(-2.7275)/(1+exp(-2.7275))
## [1] 0.06137001
 #which compares to density as..
 gden(g)
## [1] 0.06136821
#How good is the model?
 try(plot(gof(mod.rand), main=NULL))

 #add some covariates...
 mod.homoph<-ergm(g~edges+nodematch("sex")+nodematch("grade")
                          +nodematch('smokes'))  
## Evaluating log-likelihood at the estimate.
 summary(mod.homoph)
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes")
## 
## Iterations:  7 out of 20 
## 
## Monte Carlo MLE Results:
##                  Estimate Std. Error MCMC %  p-value    
## edges             -4.5787     0.1850      0  < 1e-04 ***
## nodematch.sex      0.4621     0.1311      0 0.000428 ***
## nodematch.grade    3.0493     0.1422      0  < 1e-04 ***
## nodematch.smokes   0.4062     0.1483      0 0.006196 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 6890  on 4970  degrees of freedom
##  Residual Deviance: 1712  on 4966  degrees of freedom
##  
## AIC: 1720    BIC: 1746    (Smaller is better.)
 mod.gof=gof(mod.homoph)
 plot(mod.gof, main=NULL)

  #add some simple volume bitscovariates...
 mod.p1<-ergm(g~edges+nodematch("sex")+nodematch("grade")
                  +nodematch('smokes')+sender+receiver)  
## Observed statistic(s) sender4, sender5, sender15, sender17, sender42, sender57, sender63, receiver4, receiver17, and receiver44 are at their smallest attainable values. Their coefficients will be fixed at -Inf.
## Evaluating log-likelihood at the estimate.
 summary(mod.p1)
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") + 
##     sender + receiver
## 
## Iterations:  7 out of 20 
## 
## Monte Carlo MLE Results:
##                  Estimate Std. Error MCMC % p-value    
## edges            -6.54235    1.39273      0  <1e-04 ***
## nodematch.sex     0.72340    0.15509      0  <1e-04 ***
## nodematch.grade   3.60634    0.17189      0  <1e-04 ***
## nodematch.smokes  1.27256    0.20774      0  <1e-04 ***
## sender2           2.22577    1.25583      0  0.0764 .  
## sender3           0.85358    1.31365      0  0.5159    
## sender4              -Inf    0.00000      0  <1e-04 ***
## sender5              -Inf    0.00000      0  <1e-04 ***
## sender6           1.17264    1.29221      0  0.3642    
## sender7           1.82007    1.28618      0  0.1571    
## sender8           2.33419    1.27155      0  0.0665 .  
## sender9           2.03254    1.27211      0  0.1102    
## sender10         -0.12412    1.41704      0  0.9302    
## sender11         -0.74302    1.60190      0  0.6428    
## sender12         -1.01392    1.57899      0  0.5208    
## sender13          1.44915    1.30637      0  0.2674    
## sender14          1.33064    1.32471      0  0.3152    
## sender15             -Inf    0.00000      0  <1e-04 ***
## sender16          2.34335    1.25453      0  0.0618 .  
## sender17             -Inf    0.00000      0  <1e-04 ***
## sender18          1.99117    1.29061      0  0.1229    
## sender19         -0.07322    1.40150      0  0.9583    
## sender20          1.52207    1.28439      0  0.2361    
## sender21          2.09979    1.26439      0  0.0968 .  
## sender22          0.88576    1.34555      0  0.5104    
## sender23          1.46870    1.28751      0  0.2540    
## sender24          0.88576    1.34555      0  0.5104    
## sender25          1.25861    1.34235      0  0.3485    
## sender26          1.34868    1.27622      0  0.2907    
## sender27          0.51817    1.35188      0  0.7015    
## sender28          0.30073    1.35208      0  0.8240    
## sender29          2.43662    1.28473      0  0.0579 .  
## sender30          1.36028    1.32210      0  0.3036    
## sender31          1.88595    1.28043      0  0.1408    
## sender32          0.72578    1.32209      0  0.5831    
## sender33          2.09444    1.25897      0  0.0963 .  
## sender34          2.92916    1.26585      0  0.0207 *  
## sender35          1.48498    1.28646      0  0.2484    
## sender36          2.04011    1.27129      0  0.1086    
## sender37          1.05055    1.32819      0  0.4290    
## sender38          2.06804    1.28481      0  0.1075    
## sender39          1.03670    1.40765      0  0.4615    
## sender40          1.71434    1.26586      0  0.1757    
## sender41         -0.86750    1.56951      0  0.5805    
## sender42             -Inf    0.00000      0  <1e-04 ***
## sender43          1.97972    1.28468      0  0.1234    
## sender44          1.84499    1.28992      0  0.1527    
## sender45          0.37849    1.35551      0  0.7801    
## sender46          1.20985    1.30077      0  0.3524    
## sender47         -0.07322    1.40150      0  0.9583    
## sender48          0.83024    1.30801      0  0.5256    
## sender49          1.82994    1.29996      0  0.1593    
## sender50          0.39527    1.35435      0  0.7704    
## sender51          2.59763    1.28047      0  0.0425 *  
## sender52          2.39671    1.28765      0  0.0628 .  
## sender53          0.74770    1.30858      0  0.5678    
## sender54          0.63491    1.34110      0  0.6359    
## sender55          0.87978    1.41662      0  0.5346    
## sender56          1.10414    1.30234      0  0.3966    
## sender57             -Inf    0.00000      0  <1e-04 ***
## sender58          2.32900    1.29162      0  0.0714 .  
## sender59          0.96191    1.31592      0  0.4648    
## sender60          1.20635    1.29933      0  0.3532    
## sender61          2.62317    1.27827      0  0.0402 *  
## sender62         -0.79590    1.57732      0  0.6139    
## sender63             -Inf    0.00000      0  <1e-04 ***
## sender64          1.37381    1.28712      0  0.2859    
## sender65          2.34097    1.29990      0  0.0718 .  
## sender66          2.35099    1.28906      0  0.0682 .  
## sender67          0.63782    1.41165      0  0.6514    
## sender68          0.77762    1.32296      0  0.5567    
## sender69          1.78102    1.27294      0  0.1618    
## sender70          0.09125    1.43162      0  0.9492    
## sender71          0.86733    1.45067      0  0.5499    
## receiver2        -1.43489    1.07362      0  0.1814    
## receiver3        -1.69160    1.06305      0  0.1116    
## receiver4            -Inf    0.00000      0  <1e-04 ***
## receiver5        -0.79045    0.97874      0  0.4193    
## receiver6        -1.66543    1.06255      0  0.1171    
## receiver7        -0.66723    0.95563      0  0.4851    
## receiver8         0.22141    0.91166      0  0.8081    
## receiver9         0.36172    0.86484      0  0.6758    
## receiver10       -0.87903    0.95032      0  0.3550    
## receiver11       -1.25341    1.00146      0  0.2108    
## receiver12       -2.50486    1.27755      0  0.0500 *  
## receiver13       -1.65921    1.08480      0  0.1262    
## receiver14       -0.19915    0.93421      0  0.8312    
## receiver15       -0.21006    0.90663      0  0.8168    
## receiver16       -1.03081    0.99540      0  0.3005    
## receiver17           -Inf    0.00000      0  <1e-04 ***
## receiver18       -1.00664    1.02771      0  0.3274    
## receiver19       -0.94528    0.92624      0  0.3075    
## receiver20        0.13942    0.88488      0  0.8748    
## receiver21       -1.59147    1.08360      0  0.1420    
## receiver22       -1.73458    1.09012      0  0.1116    
## receiver23       -0.12357    0.88993      0  0.8896    
## receiver24       -1.73458    1.09012      0  0.1116    
## receiver25       -0.42157    0.98180      0  0.6677    
## receiver26       -1.22260    0.96921      0  0.2072    
## receiver27        0.32993    0.87520      0  0.7062    
## receiver28       -2.47186    1.27910      0  0.0534 .  
## receiver29        0.47522    0.91342      0  0.6029    
## receiver30        0.12959    0.91004      0  0.8868    
## receiver31        0.51347    0.87020      0  0.5552    
## receiver32       -1.66313    1.06639      0  0.1189    
## receiver33        0.31993    0.87070      0  0.7133    
## receiver34        1.66643    0.85103      0  0.0503 .  
## receiver35        0.12686    0.87443      0  0.8847    
## receiver36        0.57683    0.85120      0  0.4980    
## receiver37       -0.54970    0.94518      0  0.5609    
## receiver38        0.20867    0.91396      0  0.8194    
## receiver39       -0.30300    0.98849      0  0.7592    
## receiver40       -1.12114    0.97563      0  0.2506    
## receiver41       -1.32436    0.96960      0  0.1720    
## receiver42        0.75364    0.87979      0  0.3917    
## receiver43       -0.14799    0.92850      0  0.8734    
## receiver44           -Inf    0.00000      0  <1e-04 ***
## receiver45       -1.25725    0.99393      0  0.2060    
## receiver46       -0.15188    0.90265      0  0.8664    
## receiver47       -0.94528    0.92624      0  0.3075    
## receiver48       -0.34516    0.87911      0  0.6946    
## receiver49       -0.28730    0.97578      0  0.7684    
## receiver50       -0.47568    0.91303      0  0.6024    
## receiver51        1.82942    0.84156      0  0.0298 *  
## receiver52        0.59043    0.90601      0  0.5146    
## receiver53       -1.74439    1.05338      0  0.0978 .  
## receiver54        0.41471    0.84789      0  0.6248    
## receiver55       -1.61996    1.28711      0  0.2082    
## receiver56       -1.72440    1.07485      0  0.1087    
## receiver57       -1.70193    1.08087      0  0.1154    
## receiver58        0.05836    0.94460      0  0.9507    
## receiver59       -0.63772    0.95194      0  0.5029    
## receiver60        0.34174    0.86306      0  0.6921    
## receiver61        1.19676    0.87435      0  0.1711    
## receiver62       -1.68211    1.08277      0  0.1204    
## receiver63       -0.43936    1.06709      0  0.6806    
## receiver64       -2.49890    1.28600      0  0.0521 .  
## receiver65       -0.37488    1.01075      0  0.7107    
## receiver66        0.68781    0.90364      0  0.4466    
## receiver67       -0.02045    0.94255      0  0.9827    
## receiver68       -0.77370    0.94078      0  0.4109    
## receiver69        0.15329    0.88341      0  0.8622    
## receiver70       -1.75929    1.08335      0  0.1045    
## receiver71       -0.47914    1.11343      0  0.6670    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 6890  on 4970  degrees of freedom
##  Residual Deviance: 2292  on 4826  degrees of freedom
##  
## AIC: 2580    BIC: 3517    (Smaller is better.) 
## 
##  Warning: The following terms have infinite coefficient estimates:
##   sender4 sender5 sender15 sender17 sender42 sender57 sender63 receiver4 receiver17 receiver44
 mod.gof=gof(mod.p1)
 plot(mod.gof, main=NULL)

 #model predicts better, but huge BIC cost!  Simplify?
 
 
 
 g %v% "indegree" <- indegree
 g %v% "outdegree" <- outdegree

 mod.p2<-ergm(g~edges+nodematch("sex")+nodematch("grade")
              +nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
              +mutual)
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 20:
## Optimizing with step length 0.789958896639935.
## The log-likelihood improved by 3.242.
## Iteration 2 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.1611.
## Step length converged once. Increasing MCMC sample size.
## Iteration 3 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.02584.
## Step length converged twice. Stopping.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
 summary(mod.p2)
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") + 
##     nodeicov("indegree") + nodeocov("outdegree") + mutual
## 
## Iterations:  3 out of 20 
## 
## Monte Carlo MLE Results:
##                    Estimate Std. Error MCMC % p-value    
## edges              -8.33888    0.37137      0  <1e-04 ***
## nodematch.sex       0.53239    0.12932      0  <1e-04 ***
## nodematch.grade     2.60143    0.15791      0  <1e-04 ***
## nodematch.smokes    0.86950    0.15439      0  <1e-04 ***
## nodeicov.indegree   0.29062    0.02828      0  <1e-04 ***
## nodeocov.outdegree  0.32907    0.03489      0  <1e-04 ***
## mutual              2.32249    0.23160      0  <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 6890  on 4970  degrees of freedom
##  Residual Deviance: 1321  on 4963  degrees of freedom
##  
## AIC: 1335    BIC: 1381    (Smaller is better.)
 mod.gof=gof(mod.p2)
 plot(mod.gof, main=NULL)

 mod.gendmut<-ergm(g~edges+nodematch("sex")+nodematch("grade")
              +nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
              +mutual(same="sex",diff=TRUE))
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 20:
## Optimizing with step length 0.953219017725786.
## The log-likelihood improved by 3.079.
## Iteration 2 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.3219.
## Step length converged once. Increasing MCMC sample size.
## Iteration 3 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.03451.
## Step length converged twice. Stopping.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
 summary(mod.gendmut)
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") + 
##     nodeicov("indegree") + nodeocov("outdegree") + mutual(same = "sex", 
##     diff = TRUE)
## 
## Iterations:  3 out of 20 
## 
## Monte Carlo MLE Results:
##                    Estimate Std. Error MCMC % p-value    
## edges              -8.18801    0.38824      0  <1e-04 ***
## nodematch.sex      -0.14310    0.17442      0   0.412    
## nodematch.grade     2.83441    0.15881      0  <1e-04 ***
## nodematch.smokes    0.95223    0.16382      0  <1e-04 ***
## nodeicov.indegree   0.30719    0.02848      0  <1e-04 ***
## nodeocov.outdegree  0.35012    0.03357      0  <1e-04 ***
## mutual.same.sex.1   2.77943    0.37059      0  <1e-04 ***
## mutual.same.sex.2   2.53856    0.34402      0  <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 6890  on 4970  degrees of freedom
##  Residual Deviance: 1339  on 4962  degrees of freedom
##  
## AIC: 1355    BIC: 1407    (Smaller is better.)
 #let's try a transitivity term...
# mod.tran<-ergm(g~edges+nodematch("sex")+nodematch("grade")
#                   +nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
#                                +mutual+triangle)
# summary(mod.gendmut)
 
 #turns out that won't run...lets cheat!
 mod.tran<-ergm(g~edges+nodematch("sex")+nodematch("grade")
                +nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
                +mutual+triangle, control=control.ergm(MCMLE.maxit=2))
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 2:
## Optimizing with step length 0.312988763891207.
## The log-likelihood improved by 5.029.
## Iteration 2 of at most 2:
## Optimizing with step length 0.000265821143928091.
## The log-likelihood improved by 0.2442.
## MCMLE estimation did not converge after 2 iterations. The estimated coefficients may not be accurate. Estimation may be resumed by passing the coefficients as initial values; see 'init' under ?control.ergm for details.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
 summary(mod.tran)
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") + 
##     nodeicov("indegree") + nodeocov("outdegree") + mutual + triangle
## 
## Iterations:  2 out of 2 
## 
## Monte Carlo MLE Results:
##                     Estimate Std. Error MCMC %  p-value    
## edges              -8.025142   0.405900     14  < 1e-04 ***
## nodematch.sex       0.475142   0.487653      2 0.329933    
## nodematch.grade     2.195880   0.634579      2 0.000544 ***
## nodematch.smokes    0.852256   0.583772      3 0.144378    
## nodeicov.indegree   0.211418   0.182808      1 0.247531    
## nodeocov.outdegree  0.248641   0.162040      1 0.124985    
## mutual              1.557091   0.344459     15  < 1e-04 ***
## triangle            0.106772   0.001898      7  < 1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 6890  on 4970  degrees of freedom
##  Residual Deviance: 1468  on 4962  degrees of freedom
##  
## AIC: 1484    BIC: 1536    (Smaller is better.)
 mcmc.diagnostics(mod.tran)
## Sample statistics summary:
## 
## Iterations = 16384:1063936
## Thinning interval = 1024 
## Number of chains = 1 
## Sample size per chain = 1024 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                        Mean      SD Naive SE Time-series SE
## edges                4219.3   58.02   1.8130         56.161
## nodematch.sex        2065.4   31.12   0.9726         29.214
## nodematch.grade       613.5   10.06   0.3143          8.183
## nodematch.smokes     2870.7   46.91   1.4659         44.743
## nodeicov.indegree   18093.8  128.82   4.0257        112.978
## nodeocov.outdegree  18390.6  130.86   4.0894        112.706
## mutual               2175.4   29.10   0.9094         30.775
## triangle           395589.9 7709.43 240.9196       7509.066
## 
## 2. Quantiles for each variable:
## 
##                      2.5%      25%    50%      75%  97.5%
## edges                4118   4132.5   4253   4255.0   4258
## nodematch.sex        2010   2023.2   2083   2084.0   2088
## nodematch.grade       595    606.5    619    619.2    621
## nodematch.smokes     2789   2801.2   2898   2899.0   2903
## nodeicov.indegree   17861  17919.2  18165  18173.0  18190
## nodeocov.outdegree  18155  18194.0  18463  18471.0  18490
## mutual               2126   2129.0   2193   2193.0   2194
## triangle           382591 382775.5 400279 400281.0 400299
## 
## 
## Sample statistics cross-correlations:
##                        edges nodematch.sex nodematch.grade
## edges              1.0000000     0.9992569       0.9932355
## nodematch.sex      0.9992569     1.0000000       0.9945745
## nodematch.grade    0.9932355     0.9945745       1.0000000
## nodematch.smokes   0.9998830     0.9992805       0.9935824
## nodeicov.indegree  0.9988885     0.9985409       0.9934550
## nodeocov.outdegree 0.9989202     0.9982089       0.9909539
## mutual             0.9978185     0.9959606       0.9883397
## triangle           0.9974954     0.9953521       0.9866779
##                    nodematch.smokes nodeicov.indegree nodeocov.outdegree
## edges                     0.9998830         0.9988885          0.9989202
## nodematch.sex             0.9992805         0.9985409          0.9982089
## nodematch.grade           0.9935824         0.9934550          0.9909539
## nodematch.smokes          1.0000000         0.9987943          0.9987190
## nodeicov.indegree         0.9987943         1.0000000          0.9975310
## nodeocov.outdegree        0.9987190         0.9975310          1.0000000
## mutual                    0.9971929         0.9950552          0.9966797
## triangle                  0.9968425         0.9944679          0.9963066
##                       mutual  triangle
## edges              0.9978185 0.9974954
## nodematch.sex      0.9959606 0.9953521
## nodematch.grade    0.9883397 0.9866779
## nodematch.smokes   0.9971929 0.9968425
## nodeicov.indegree  0.9950552 0.9944679
## nodeocov.outdegree 0.9966797 0.9963066
## mutual             1.0000000 0.9998440
## triangle           0.9998440 1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##              edges nodematch.sex nodematch.grade nodematch.smokes
## Lag 0    1.0000000     1.0000000       1.0000000        1.0000000
## Lag 1024 0.9981101     0.9977836       0.9970511        0.9980660
## Lag 2048 0.9960372     0.9954138       0.9941653        0.9959323
## Lag 3072 0.9938211     0.9929642       0.9914632        0.9936546
## Lag 4096 0.9914587     0.9904576       0.9886689        0.9912237
## Lag 5120 0.9889850     0.9879250       0.9857973        0.9886773
##          nodeicov.indegree nodeocov.outdegree    mutual  triangle
## Lag 0            1.0000000          1.0000000 1.0000000 1.0000000
## Lag 1024         0.9974614          0.9973679 0.9982533 0.9983115
## Lag 2048         0.9948052          0.9948389 0.9963786 0.9964605
## Lag 3072         0.9921616          0.9923010 0.9943711 0.9944328
## Lag 4096         0.9894583          0.9896352 0.9922170 0.9922461
## Lag 5120         0.9866739          0.9869925 0.9899313 0.9899126
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##              edges      nodematch.sex    nodematch.grade 
##            -185.35            -112.50             -33.57 
##   nodematch.smokes  nodeicov.indegree nodeocov.outdegree 
##            -245.53            -117.93             -99.36 
##             mutual           triangle 
##            -279.57           -7890.08 
## 
## Individual P-values (lower = worse):
##              edges      nodematch.sex    nodematch.grade 
##       0.000000e+00       0.000000e+00      4.930723e-247 
##   nodematch.smokes  nodeicov.indegree nodeocov.outdegree 
##       0.000000e+00       0.000000e+00       0.000000e+00 
##             mutual           triangle 
##       0.000000e+00       0.000000e+00 
## Joint P-value (lower = worse):  1.491166e-87 .
## Package latticeExtra is not installed. Falling back on coda's default MCMC diagnostic plots.

## 
## MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).
 test=simulate.ergm(mod.tran,1)
 plot(test)

 mod.gwesp<-ergm(g~edges+nodematch("sex")+nodematch("grade")
                +nodematch('smokes')+nodeicov("indegree")+nodeocov("outdegree")
                +mutual+gwesp)
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 20:
## Optimizing with step length 0.564907830439016.
## The log-likelihood improved by 3.669.
## Iteration 2 of at most 20:
## Optimizing with step length 0.952429948637443.
## The log-likelihood improved by 2.673.
## Iteration 3 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.1513.
## Step length converged once. Increasing MCMC sample size.
## Iteration 4 of at most 20:
## Optimizing with step length 1.
## The log-likelihood improved by 0.05256.
## Step length converged twice. Stopping.
## Evaluating log-likelihood at the estimate. Using 20 bridges: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .
## This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
 summary(mod.gwesp)
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ edges + nodematch("sex") + nodematch("grade") + nodematch("smokes") + 
##     nodeicov("indegree") + nodeocov("outdegree") + mutual + gwesp
## 
## Iterations:  4 out of 20 
## 
## Monte Carlo MLE Results:
##                    Estimate Std. Error MCMC % p-value    
## edges              -8.02670    0.39054      0  <1e-04 ***
## nodematch.sex       0.50762    0.11943      0  <1e-04 ***
## nodematch.grade     2.20383    0.17854      0  <1e-04 ***
## nodematch.smokes    0.78078    0.14779      0  <1e-04 ***
## nodeicov.indegree   0.23866    0.02982      0  <1e-04 ***
## nodeocov.outdegree  0.27343    0.03474      0  <1e-04 ***
## mutual              2.05411    0.24294      0  <1e-04 ***
## gwesp               0.64075    0.14652      0  <1e-04 ***
## gwesp.decay         0.06569    0.15439      0   0.671    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##      Null Deviance: 6890  on 4970  degrees of freedom
##  Residual Deviance: 1297  on 4961  degrees of freedom
##  
## AIC: 1315    BIC: 1373    (Smaller is better.)
 mcmc.diagnostics(mod.gwesp)
## Sample statistics summary:
## 
## Iterations = 16384:4209664
## Thinning interval = 1024 
## Number of chains = 1 
## Sample size per chain = 4096 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                        Mean      SD Naive SE Time-series SE
## edges               -2.1826  18.993   0.2968         1.0084
## nodematch.sex       -1.2117  14.086   0.2201         0.7904
## nodematch.grade     -2.1377  14.659   0.2290         0.8521
## nodematch.smokes    -0.9524  15.984   0.2497         0.8639
## nodeicov.indegree   -9.9270 110.148   1.7211         5.7165
## nodeocov.outdegree -19.4607 108.374   1.6933         5.6268
## mutual              -1.3481   7.733   0.1208         0.4565
## gwesp               -3.3380  22.683   0.3544         1.2833
## gwesp.decay         -1.6761  13.610   0.2127         0.7583
## 
## 2. Quantiles for each variable:
## 
##                       2.5%    25%     50%    75%  97.5%
## edges               -39.00 -15.00  -2.000 11.000  35.00
## nodematch.sex       -29.00 -11.00  -1.000  9.000  26.00
## nodematch.grade     -30.00 -12.00  -2.000  8.000  27.00
## nodematch.smokes    -31.00 -12.00  -1.000 10.000  31.00
## nodeicov.indegree  -225.00 -87.00  -9.000 65.000 210.00
## nodeocov.outdegree -230.62 -91.00 -20.000 54.000 197.00
## mutual              -17.00  -6.00  -1.000  4.000  14.00
## gwesp               -47.02 -18.36  -3.720 11.907  43.07
## gwesp.decay         -27.20 -10.88  -2.058  7.213  26.51
## 
## 
## Sample statistics cross-correlations:
##                        edges nodematch.sex nodematch.grade
## edges              1.0000000     0.7920421       0.8285412
## nodematch.sex      0.7920421     1.0000000       0.6338130
## nodematch.grade    0.8285412     0.6338130       1.0000000
## nodematch.smokes   0.8821243     0.6993325       0.7473641
## nodeicov.indegree  0.9243292     0.7224975       0.6880821
## nodeocov.outdegree 0.9485801     0.7264189       0.7353055
## mutual             0.7850074     0.6248356       0.7928949
## gwesp              0.9225984     0.7214260       0.8427424
## gwesp.decay        0.7798153     0.5812022       0.7877265
##                    nodematch.smokes nodeicov.indegree nodeocov.outdegree
## edges                     0.8821243         0.9243292          0.9485801
## nodematch.sex             0.6993325         0.7224975          0.7264189
## nodematch.grade           0.7473641         0.6880821          0.7353055
## nodematch.smokes          1.0000000         0.7572631          0.8029944
## nodeicov.indegree         0.7572631         1.0000000          0.8845197
## nodeocov.outdegree        0.8029944         0.8845197          1.0000000
## mutual                    0.7077234         0.7101292          0.7345627
## gwesp                     0.8103180         0.8669983          0.8927071
## gwesp.decay               0.6894998         0.7446771          0.7637741
##                       mutual     gwesp gwesp.decay
## edges              0.7850074 0.9225984   0.7798153
## nodematch.sex      0.6248356 0.7214260   0.5812022
## nodematch.grade    0.7928949 0.8427424   0.7877265
## nodematch.smokes   0.7077234 0.8103180   0.6894998
## nodeicov.indegree  0.7101292 0.8669983   0.7446771
## nodeocov.outdegree 0.7345627 0.8927071   0.7637741
## mutual             1.0000000 0.8189673   0.7645249
## gwesp              0.8189673 1.0000000   0.8506887
## gwesp.decay        0.7645249 0.8506887   1.0000000
## 
## Sample statistics auto-correlation:
## Chain 1 
##              edges nodematch.sex nodematch.grade nodematch.smokes
## Lag 0    1.0000000     1.0000000       1.0000000        1.0000000
## Lag 1024 0.7356949     0.7243649       0.8194960        0.7344987
## Lag 2048 0.5964351     0.5895078       0.7048158        0.6006431
## Lag 3072 0.5078308     0.4968075       0.6115425        0.5079153
## Lag 4096 0.4374859     0.4333647       0.5375935        0.4377303
## Lag 5120 0.3812865     0.3798584       0.4708243        0.3895060
##          nodeicov.indegree nodeocov.outdegree    mutual     gwesp
## Lag 0            1.0000000          1.0000000 1.0000000 1.0000000
## Lag 1024         0.7273865          0.7310531 0.8378584 0.7870712
## Lag 2048         0.5831148          0.5891169 0.7203466 0.6655837
## Lag 3072         0.4913930          0.4989458 0.6233013 0.5797273
## Lag 4096         0.4242818          0.4363518 0.5490370 0.5117562
## Lag 5120         0.3672405          0.3829011 0.4868593 0.4542104
##          gwesp.decay
## Lag 0      1.0000000
## Lag 1024   0.7811626
## Lag 2048   0.6533085
## Lag 3072   0.5649845
## Lag 4096   0.4998976
## Lag 5120   0.4413901
## 
## Sample statistics burn-in diagnostic (Geweke):
## Chain 1 
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##              edges      nodematch.sex    nodematch.grade 
##             0.5028             1.3126             1.1333 
##   nodematch.smokes  nodeicov.indegree nodeocov.outdegree 
##             0.2698             0.2840             0.3324 
##             mutual              gwesp        gwesp.decay 
##             1.5353             0.8789             0.8636 
## 
## Individual P-values (lower = worse):
##              edges      nodematch.sex    nodematch.grade 
##          0.6151170          0.1893017          0.2570862 
##   nodematch.smokes  nodeicov.indegree nodeocov.outdegree 
##          0.7872953          0.7764182          0.7396138 
##             mutual              gwesp        gwesp.decay 
##          0.1246995          0.3794302          0.3877901 
## Joint P-value (lower = worse):  0.09919254 .
## Package latticeExtra is not installed. Falling back on coda's default MCMC diagnostic plots.

## 
## MCMC diagnostics shown here are from the last round of simulation, prior to computation of final parameter estimates. Because the final estimates are refinements of those used for this simulation run, these diagnostics may understate model performance. To directly assess the performance of the final model on in-model statistics, please use the GOF command: gof(ergmFitObject, GOF=~model).
 test=simulate.ergm(mod.gwesp,1)
 plot(test)

 try(plot(gof(mod.gwesp), main=NULL))

 ################################################################
 #latent space models.  These models allow us to fit a network without 
 #having to specify the actual tie formation patterns...
 #code cribbed from Jake Fisher
 #################################################################
 
 # The most user-friendly latent space model software is in the latentnet package
 # (more recent models are provided by the amen package).
 #install.packages("latentnet")
 
 library(latentnet)
## 
## latentnet: version 2.8.0, created on 2017-08-20
## Copyright (c) 2017, Pavel N. Krivitsky, University of Wollongong
##                     Mark S. Handcock, University of California -- Los Angeles
##                     with contributions from
##                     Susan M. Shortreed
##                     Jeremy Tantrum
##                     Peter D. Hoff
##                     Li Wang
##                     Kirk Li, University of Washington
##                     Jake Fisher
##                     Jordan T. Bates
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("latentnet").
## NOTE: BIC calculation prior to latentnet 2.7.0 had a bug in the calculation of the effective number of parameters. See help(summary.ergmm) for details.
## NOTE: Prior to version 2.8.0, handling of fixed effects for directed networks had a bug. See help("ergmm-terms") for details.
start.time <- Sys.time()
  latent.fit <- ergmm(g ~ euclidean(d = 2))
runtime=Sys.time()-start.time;
runtime;
## Time difference of 1.408969 mins
 summary(latent.fit)
## NOTE: It is not certain whether it is appropriate to use latentnet's BIC to select latent space dimension, whether or not to include actor-specific random effects, and to compare clustered models with the unclustered model.
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ euclidean(d = 2)
## Attribute: edges
## Model:     Bernoulli 
## MCMC sample of size 4000, draws are 10 iterations apart, after burnin of 10000 iterations.
## Covariate coefficients posterior means:
##             Estimate   2.5%  97.5% 2*min(Pr(>0),Pr(<0))    
## (Intercept)   1.7690 1.3795 2.1624            < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Overall BIC:        1910.635 
## Likelihood BIC:     1127.11 
## Latent space/clustering BIC:     783.5258 
## 
## Covariate coefficients MKL:
##              Estimate
## (Intercept) 0.9692998
 plot(latent.fit)

 mcmc.diagnostics(latent.fit)
## Chain 1 
## Lag 0 
##              lpY    beta.1     Z.1.1     Z.1.2
## lpY    1.0000000 0.4594187 0.1319291 0.1317427
## beta.1 0.4594187 1.0000000 0.1526472 0.2620373
## Z.1.1  0.1319291 0.1526472 1.0000000 0.2264011
## Z.1.2  0.1317427 0.2620373 0.2264011 1.0000000
## 
## Lag 10 
##               lpY     beta.1      Z.1.1      Z.1.2
## lpY    0.45622566 0.29748463 0.10352657 0.07067086
## beta.1 0.30626407 0.31471132 0.06964949 0.06954152
## Z.1.1  0.09422899 0.05918816 0.51834688 0.09056140
## Z.1.2  0.07853141 0.06887424 0.10328924 0.45935555

## [[1]]
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.0125
## Probability (s) = 0.95 
##                                               
##         Burn-in  Total Lower bound  Dependence
##         (M)      (N)   (Nmin)       factor (I)
##  lpY    40       7800  600          13.0      
##  beta.1 30       7170  600          12.0      
##  Z.1.1  40       7970  600          13.3      
##  Z.1.2  100      17340 600          28.9
 plot(gof(latent.fit), main=NULL)

 # Can add additional dimensions...
 start.time <- Sys.time()
  latent.fit <- ergmm(g ~ euclidean(d = 3))
runtime=Sys.time()-start.time;
runtime;
## Time difference of 1.981795 mins
plot(gof(latent.fit), main=NULL)

 # ... latent groups ...
start.time <- Sys.time()
 latent.fit <- ergmm(g ~ euclidean(d = 2, G = 5))
runtime=Sys.time()-start.time;
runtime;
## Time difference of 2.776853 mins
plot(latent.fit)

plot(gof(latent.fit), main=NULL)

 # ... or homophily effects
start.time <- Sys.time()
 latent.fit <- ergmm(g ~ nodematch("grade") + nodematch("sex")+euclidean(d = 2))
runtime=Sys.time()-start.time;
runtime
## Time difference of 1.369521 mins
summary(latent.fit)
## 
## ==========================
## Summary of model fit
## ==========================
## 
## Formula:   g ~ nodematch("grade") + nodematch("sex") + euclidean(d = 2)
## Attribute: edges
## Model:     Bernoulli 
## MCMC sample of size 4000, draws are 10 iterations apart, after burnin of 10000 iterations.
## Covariate coefficients posterior means:
##                 Estimate     2.5%   97.5% 2*min(Pr(>0),Pr(<0))    
## (Intercept)     -1.11165 -1.58619 -0.6259               <2e-16 ***
## nodematch.grade  3.67247  3.23129  4.0931               <2e-16 ***
## nodematch.sex    0.27219 -0.20578  0.6836                0.223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Overall BIC:        1741.307 
## Likelihood BIC:     1161.369 
## Latent space/clustering BIC:     579.9382 
## 
## Covariate coefficients MKL:
##                   Estimate
## (Intercept)     -1.9357062
## nodematch.grade  3.2004359
## nodematch.sex    0.2859632